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1 Objectives 


A well-documented problem in the analysis of data collected WjP^ cecr ^ 

r r » 

data dropouts an irregularly-spaced sampling grid, and time-delayed sam 
ptog " d“a insularities render most traditional si^ra > P ™g 

techniques unusable, and thus the data must be interpo ^ 

- , J fnrp scientific analysis techniques can be applied. In addition, tne ex 
lie Of data collected by scientific instrumentation presents 
many challenging problems in the area of compression, visualization, and 
ai^vsis 'Therefore, a representation of the data is needed which provides 
rllure which is conducive to these applications. Wavelet represen a 
tions of data have already been shown to possess excellent cl.aractenstics 

“TrlltlrJK “develop a new adaptive Altering al- 
gorithm for image restoration and compression. The algont ™ b |‘j 

tow computational complexity and a fast convergence rate. Th^ will make 
the algorithm suitable for real-time applications. The algonthm should 
able to remove additive noise and reconstruct lost data samples rom images. 

2 Introduction & Accomplishments 

The problems of image restoration and compression haw revived witepre^ 
attention in the past two decades. There are several linear and nonlinea 
technioues that have been proposed. One of the methods used is linear pre- 
diction using adaptive filtering. The goal of the adapt. vefilter *to minimize 
a cost function in real time with low numerical complexity and simplicity 
of implementation. Most of the adaptive algorithms have to trade-off b 
tween computational complexity and performance. For instance, the leas 
mean sauared algorithm (LMS) is widely used in many applications becau 
of its tow computationally cost even though it has a somewhat dow ^nw- 
eence rate In contrast, the recursive least square algorithm (RLS) has las 
Tn^nce and tow mean squared error, but it has a high computational 

COS Recentlv there has been some interest in the implementation of adap- 
tive filter! in the frequency domain. Advantages in performance have been 
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reported in several applications. In this project we ^develop ns *^aptive ah 

construction for use in image restoration and compression. I P 
a r-nrtm nl i sh cd the following in this project. 


2.1 Accomplishments 

1 Developed a new adaptive filtering algorithm called the Fast Euclidean 
' Direction Search (FEDS) method. This algorithm has a tow 
tional complexity of 0(N ) and a convergence rate which comp 

to the RLS. 

2. Developed 2-D versions of the FEDS algorithm and their variations, 

3 Investigated the use of the 1-D and 2-D FEDS algorithms for many 
' applications. The results are quite impressive and outperform thos 

existing methods. 

4. Implemented FEDS algorithms in the wavelet domain, 
r implemented an algorithm for lost-sample-recovery in conjunction with 

' the FEDS method. We used this for image restoration and compression 

with encouraging results. 

Extra Work We have started working on the design and implementation of 
2 D multiplierless filters. Multiplierless filters are those that have coe cien 
wM c ra" P owers-of-two. Therefore, multipliers are simply replaced by 
vff „_ These filters have high computational speed and low cost. Multi 
shdt rs. for use in image restoration, 

£££ and videTp^ceU. vl have developed some design methods 
for multiplierless filters. Further design and implementation work 
topic will be conducted in the future. 


2.2 Publications 


The following publications resulted from this project. 
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The rest of the report is organized as follows. In Section 3, the deriva- 
tion of the FEDS algorithm is presented. In Section 4, the wavelet domain 
implementation of the FEDS is described. Section 5 reports some results 
the application of these algorithms for image restoration. In Section 6, some 


4 


new results are reported on the design of multiplierless digital filters. Section 
7 is the conclusion. 


3 Euclidean Direction Search Adaptive Algo- 
rithm 

The Euclidean Direction Search (EDS) algorithm has been reported by the 
PI and his research team in [6], [7], [8], [9], [10]. In this section, we present 

the basic derivation and some features of the algorithm. 

The EDS algorithm is essentially a derivative of the Powell and Zangwill 
DS algorithm. In adaptive applications, we consider the following quadratic 

optimization problem, 


J n (w) = min {w T Q (n) w - 2w T r (n) + o 1 (n)} 
v ' w eR N 


( 1 ) 


where Q (n) is an N x N symmetric positive definite matrix r (n) is a column 
vector of order N, w is the unknown weight vector, and a is a variance o^ 
the desired signal, d. The minimum of J n (w) can be calculated by setting 
the gradient to zero,V 4 ^ = 0 . Obviously, the optimal solution is 

W< %brthe least squares optimization problem, the cost function can be ex- 
pressed as, 


A n-i e 2 (n) 

i=l 

(2) 

(d(i) - w T x (i)) 2 

i 

(3) 

2=1 

a\ (n) - 2w r (n) r (n) + 


w r (n) Q (n) w (n) 

(4) 
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where 


n 

Q(n) = ^A n_i x(i)x T (i) 

i= 1 

= AQ (n - 1 ) + x(n)x r (n) 

n 

r ( n ) = 'y j ^ n ^ x ^ 

i= 1 

= Ar(n-l) + d(n)x(n). 


(5) 

( 6 ) 

(7) 

( 8 ) 


* In) = Ixln) x (» - 1) x (n - 2) ... x (n - N + l)fis an input vector to the 
adaptive filter of length JV at time n, d(n) is a desired signal, and A is a 
forgetting factor. Consider the update weight iterative function 


w (n + 1) = w (n) + a (n) g (n) 


(9) 


where a (n) and g (n) are step size and search direction at iteration n, ^re- 
spectively. One can obtain the best step size to minimize the next step 
function by setting V Q = 55 u a ' s touo ’ 

dJn (w-rog) _ q _ 2g T Q (n) (w+ag) - 2g T r ( n ) . 

r)rv 


As a result, 


a = — ■ 


g T (Q (n) w - r (n)) 


( 10 ) 


g T Q («) 8 

The proposed algorithm utilizes the Euclidean direction, that is g« = 
r n 0 10 0l T where 1 is located in the i-th position, 

rectional search allows us to update the weight vector in a simple »• 
Considering equation (10), when the i - th position of the duechon >s se- 
lected we pick up the i - th element of vector (Q (n) w - r (n)), and simp y 
divide it by Q The EDS algorithm needs totally N rounds to com- 

pletely update all weight elements in a filter. The algorithm is described 
explicitly in Table 1. We denote (w), and (r) 4 as the * - t/ielement o wag ^ 
vector, w, and, cross correlation vector, r, respectively, q i 
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matrix Q. 


For n = 1,2,... 

For i = 1 , ■ • • ) N 

(1) £ = (q (l) ) T w - (r)i ; 

(2) a=(q«) i; 

(3) if a > 0 a = - s ~ ; (w)* = ( w )t + a] 

(4) qW = Aq (l) + (x {k + 1))* x ( fc + !) 

(5) (r)» = * ( r )i + ( x + *))i d ( k + l) 

end i 


Table 1. EDS Algorithm 

The EDS algorithm needs N (3 N + 3) multiplications for each incoming 
d ajsample. In other words, it has 0(«*) 

This is relatively high compared to other algorithms, although EDb has com 
perfon^l There is however, an alternative approach which can 
reduce the computational cost to O(N). The new algorithm is the so cal e 

EDS1 or Fast-EDS (FEDS) [6], [7], [81, [9], [10]. In ° rder 
nutations to O (N), we need to modify the objective function . 
exponentially weight least square form as follows: 


fc-1 


N 


min J n (w) = + j 

w eR N i=0 j= 1 

\2l 


w ( n) T x(xiV + j)) } + 


^2 (d (kN + j) - w (n) T x (kN + j )) 
i =1 


(») 


with n = kN + 1 where k is the number of full blocks within N and ! is the 
number of samples remaining. The new cost function leads to modifications 


7 



of Q(n), r (n), and a 2 (n) as 

Q (n) = Q (^N + l) 

k-\ / N 

= Y xk ~ { [Y,*( iN+j)x ( iN +j) ) + 


i=0 


0=1 


y; x (kN + j) x (kN + j) T ; 

3= l 

r (n) = r (kN + l ) 


fc-i 


N 


= ^d(iiV + j)x(i^ + j) 1 + 


i=0 


o =1 


d, (kN + j) x (kN + j) 
j =1 

cr 2 (n) = o 2 (kN + l) 


k-\ 


N 


= E d2(zAr+j) + 


2 = 0 


0 = 1 


]Td 2 (fciV + j) 


j=i 


( 12 ) 


(13) 


(14) 


To compute 
written as, 


these parameters recursively, equation (12) and (13) can be 


Q (kN + l) = Q (n) 

= AQ (kN) + Q (kN + 0 , 

r (kN + l) = r(n) 

= Xr{kN) + r(kN + 1), 


(15) 

(16) 


where 

1 ■ T 

Q{kN + l) = E x (ZciV + j) x (kN + j) 

j = 1 

l 

r (kN + l) = ^d(fciV + j)x(fciV+j). 

J = 1 


(17) 

(18) 
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The FEDS algorithm differs from the EDS algorithm in the following ways. 
First it updates the correlation matrix, Q, and cross correlation vector r 
££ he above equations. Secondly, only one Euclidean direction search is 
conducted per incoming sample of data. The FEDS algorithm is illustrated 

in Table 2. 


For k = 1,2, ... 

Step (I): For i = 1,- • • ,N 

(1) = Aqf } ; 

( 2 ) Mi =j* (r)i ; 

(3) Q = Q + x (kN + i) x T (kN + i ) ; 

(4) r = r + d {kN + i)x {kN + i) ; 

(5) £ = (q (l) + q (i) ) T w - (r + r)i ; 

(6) a = (q^)* + (q (i) )i ; 

(7) If a > 0 Q = -f; (w) i = (w) i + a; 
end i 

Step (II) _ _ 

(1) Q = Q + Qi Q = °> 

(2) r = r + r; r = 0; 


Table 2. FEDS Algorithm 

It is shown in [7] that the computational count of the FEDS algorithm is 
4JV + 2 multiplications, that is O {N). 


3.1 Stability of The EDS Algorithm 

Consider a time-invariant quadratic optimization problem, 

j (w) = min {w r Qw — 2w T r + a 2 } ( 19 ) 

where Q is an JV x N real symmetric positive definite matrix r and variable 
w are vectors of order JV, and o is a scaler. The notation (,f means the 

‘“By' Introducing a constant M to regulate the convergence rate, the New- 
ton’s method can be expressed as: 

w (k + 1) = w(fc) - ^Q _1 V(fc), 


where w(fc) denotes the variable vector at step k, V(fc) V |w_w(fc) 
2Qw (fe) - 2r, and 0 < p. < 2. 
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The steepest descent method avoids dealing with matrix Q ' 1 and updates 
the variable vector w on the negative gradient direction as follows. 

w (k + 1) = w(fc) — — V(A:), 

where 0 < a < and A max is the largest eigenvalue of Q. 

The Euclidean direction set (EDS) method in a constant environment 
equivalent to Gauss-Seidel method, which can be written as: 

( 20 ) 


w (k + 1) = w (k) - -A x V(fc) , 


where A is an N x N triangular matrix or block triangular matrix of Q as 

shown m [17] ^ of the EDS algorithm can be obtained by defining 

an error vector as C(fc) A w(fc) - w., where w = Q-'r » the op timaj 
solution of (19). Subtracting w* from both sides of (20), and using ( ) 
2(Qw(fc) - QwJ, the error recursive equation becomes 

C (Jt +!) = (/ — A-'Q)C(k), < 21 > 

where I is the identity matrix. With an initial value of C(0), for the error 
vector, the solution of (21) is 

( 22 ) 


C(fc) = (/ - A _1 Q) fe C(0). 


-A- J Q) < 1, 


It is well known that C{k) -> 0 as k -» oo, if and only p{I 
where p(.) denotes the maximum absolute eigenvalue of (.). 

Based on the fact that matrix A is a lower triangular submatrix o Q, 
the following theorem provides a sufficient condition for the convergence o 
the EDS algorithm. It is worth noting that this sufficient condition is inde- 
pendent of the largest eigenvalue of matrix Q. 

Theorem 3.1 Let Q be an N x N symmetric positive definite matrix a 
matrix A be the lower triangular submatrix or a block triangular subma^ 
trix of Q. The euclidean direction search (EDS) algorithm w(fc + ) 

_ lA _1 V(/c) converges to the optimal solution. 

Proof. Yet B = A + A T -Q. It is easy to verify that matrix B is a block 
diagonal submatrix of Q, with each block size being 1 x 1 or 2 x % As prow 
in the next subsection, matrix B is symmetric and positive definite. Now, 
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assume that A and x are the eigenvalue and associated eigenvector of matrix 
A' 1 Q where 

(I- A _1 Q)x= (l-A)x. 

Multiplying by x*A. taking the absolute value of both sides, and rear- 
ranging gives _ |x*Ax — x*Qx| 

^ ^ |x*Ax| 

where (.)* denotes the conjugate transpose of (.)• , _ 

Recall that for x 6 C";x / 0, x'Qx > 0 , x‘Bx > 0, and Re{x Ax) - 

x*(^£)x = x*(^)x>0, 

Re{x*Ax — x*Qx} < Re{x Ax}; 

Im{x* Ax — x*Qx} = Im{x*Ax). 

That is, |1 - A| < 1, which implies that p{ I - A _1 Q) < 1- That is, C{k) -► 0 
as jfe -> oo in (22), and the conclusion follows. (j 

3.2 Convergence Rate of The EDS Method 

In the previous section, we have shown that |1 - A| < 1, which implies 
that A + 0. Therefore, using similarity transformation, we may express (22) 


as 


C(fc) = S(I - A) fc S" 1 C(0), 


where matrix A is a diagonal matrix containing the eigenvalues of matrix 
(A -1 Q) and A -1 Q = SAS' 1 . For convenience, define a new vector V(fc) 

S^Cffc), so that , 

K) V(fc) = (I — A) fc V(0). 

Clearly, the convergence rate of each element t*(fc) in vector V(fc) is 

dependent on the associated eigenvalue A, of matrix (A <4)- that 

The quantity r, = 1 - A, is known as the “geometric ratio . Note that 
when the absolute value of r, is less than 1, the rate of convergence increases 

95 ^s^wdl known, the overall convergence rate r cannot be expressed in 
a simple closed form. But fortunately, the absolute value of geometric ratio 

r is lower bounded. So, we indicate that the convergence performance of he 

EDS method is superior to the steepest descent method by showing that the 
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lower bound \r\^ nd in the EDS method is lower than that of the steepest 

d " e ho°ut the rest of the chapter, A.(.) A (-) and demote 

the ith, the largest and the smallest eigenvalues of matnx “dTr L - 
The ith geometric ratio, for the steepest descent method, is r t 
M < P < The overall convergence rate rs lower 

bounded as ^ ^ - M Q)\) ■ < 23 > 

The best step size p for the convergence occurs at the cross-point: AJQ)' 
1 = 1 - pA ml „(Q),wWch gives p = S " mg “ mt ° <23) ’ 

yields |r| > S °’ in the steepeSt dCSCent meth ° d ’ 

_ A m a x(Q) ~ ^niin(Q) 

I ^ I bound 


Amax(Q) d* A m in(Q) 
2 

1 - 


x(Q) 


(24) 


n CQ) 


+ 1 


where Amax ,^ equals the condition number of Q. 

Amin Q , pv atirr Anm AQ) decreases Hence, the 

Note that \r\ bound decreases as the ratio Amin(Q) de 

convergence rate increases as the ratio fegi decreases. 

In order to compare the bounds between two methods, let s prove the 

following lemma first. 

Lemma 3.1 Assume Q is on NxN symmetric positive definite matrix, 
and matrix B is a block diagonal submatnx of Q, such as 


9n 912 

921 922 

0 


0 


0 


933 934 

943 944 


0 


0 


0 


<7n-i,jv-i 9n,n- i 
9n,n- i 9n,n 


unth each block being square. Then, matrix B is symmetric, positive 
definite, and 


A 


r min (B) > A min (Q); A max (B) < A max (Q)- 
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Proof For * = 1, 2, • • • , N, the ith principle subrnatrix Q» of Q is the i x » 
submatrix consisting of the intersection of the first i rows and columns o^. 
let B denote the jth block submatrix m B, i.e. B - Li- 

Shrce any principle subrnatrix of a symmetric positive definite matnx .s sy.n- 

metri expositive definite and MQ) < MQ.) * WQ) ^ 
it such that 1 < k < i [14], it therefore, follows that the first b 
symmetric positive definite, and A mi „ BO > A min (Q); A max (BO < A max iWh 
y S tL there exists the permutation matrices which can br.ng the 

other blocks being the first block without effect of the symmetnc and pos.t.ve 
definite properties [16], and 


‘min 


miax 


(B) = mill {A m in (By)} ; 

(B) = max {A max (By)} • 

j=l, "L 


We thus proved that matrix B is symmetric, positive definite, and 

A min (B) > A min (Q); Ajnax(B) < A nia x(Q)- # 

In the EDS method, the overall convergence rate |r| > " Xi<yA 

and I x *Qx| 


I ^ I bound 


sup { 

A -1 Qx=Ax;x/0 


x*Ax 


}. 


( 26 ) 


Since x - Qx is real, the supremum occurs when x Ax ls “||^0Oere' 
if x- Ax is real, then x' Ax = x* ^x and x* Ax - x Qx = x -^x. There 

fore, 


I ^ I bound — 


r jx*Qx - x*Bx| 


1 - 

x€R n ;x^o 
1 - 2 
xGR N ;x^O 
2 


A 5 ^0 x*Qx + x’Bx 

x*Ax£R ^ 

2__ if x*Qx > x*Bx, 

sup 

■^ N ;x^o 

if X-Qx < x Bx; 

sup J 


1 


if x*Qx > x*Bx; 


< 


■Amax (Q) _|_ ^ 

1 - [T.L -7 i f **Q X < x * Dx: 


( 27 ) 


A min 


+i 


That is, | r 


bound 


< max < 1 — 


^el ’ 1 J' 
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Recall Lemma 3.1, and note that when the equality in (25) occurs, \r\ bomd 

(Q) an( [ A max(Q) > true j n general, which 


^ mi C ^max(Q) \ ^naax 

0* Therefore, \ ^ A m ; n iBi ^mintcj) '’mnnwi , , 

asserts that |r|2L in the EDS is lower than that in steepest descent method. 
If \naxjQj ig close to i 5 with proper step size /x, both EDS and steepest 

descent method converge fast. However, when is large, the last ‘ < 

” in (27) is a very conservative step. Since matrix B is a bloc jk diagonal 
submatrix from Q and matrix Q is central majorized, ^ **Q*+** Bx 

is most likely close to the ratio £$££$ ■—* <* " 


Amax(B) — A m ln(Q) | p wlieil 
Amax (B)+A m in(Q) ’ ' A 


X(Q) 


min(Q) 


is large, 




sup 

xeC N ;x?tO 

A m in(B) 


f |x*Ax 

l R 


Ar 


|x*Ax| 

i(Q) 


x*Qx[| 


A 


min\ 


- , ( 28 ) 

n(B) + A m in(Q) 

This implies that the convergence rate of EDS is limited by the condition 
number of matrix Q somewhat too. In addition, (28) also exp 01 s e ac 
that the double direction search seems to converge more rapidly than the 
one direction search because A min (B) is smaller and more close to A min (q) in 

EDS2 

Example 3.1: Consider a simple example of a single-input adaptive 
linear combiner with two weights. The input and desired signals are sampled 
sinusoids at the same frequency, with M = 6 samples per cycle. e ■ mpu 

correlation matrix Q and the correlation vector r are calculated as follows. 

i.5 0.25 

.25 0.5 

n T 

r = E{[d k x k 4 xfc-i]} T =[o • 

Solving det(AI - Q) = 0, gives that A min (Q) = 0.25, A max (Q) = 0.75. From 
(24), \r\ bound 1 

In EDS, A = 


Q 


E 


x% x k x k _ 

X k —iX k Xf. 


* fc -i n r°-' 

1-1 Jj L°- 5 


' 0.5 
0.25 

0 1 
0.5 j 

; B = 

O O 
bi 

0 

0.5 

= 0.75, 

Aniax 

(A- l Q) 

= 1, 

and 


= 4 based on (26). 

Amin l XV Wl — U. I u, t'maxv"- "w 1 ,utmn . i • • 11 

FVom (28) and A min (B) = A„,«(B) = 0.5, we conclude that |r| > 5 minimally. 
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4 Wavelet Transform 

In this section, wavelet transform is briefly overviewed. The simplest ap- 
proach to start is to consider one-dimensional wavelet transform. ~ 
transform is a signal decomposition method using ortho 8 on » ^ 8 

with a signal x (t), which can be expressed as a linear combination as, 

(29) 


X 


(t)=£<W, !<) 


where I is a integer index, d, are decomposition coefficients, i md * denote 
the expansion set. For wavelet analysis, the expansion set ,s constructed 
from a mother wavelet, The scaling and translating version of the mot 
wavelet, produces more orthonormal basis, which are, 


(*) = 24 (Wt - k) 

where j, k are integers. Hence equation (29) becomes 

: ( t j = y y ‘ dji.V' 1 ,,! (f) 

j k 


X I 


(30) 

(31) 

(32) 


In multiresolution analysis, there is another basic function called ‘ “f 1 ” 8 

tion V (t), in addition to wavelet i/i (t). These two functions are h.ghly related 

and’used to form subspaces in L 2 (R). The relationship of subspaces can be 
written as, ~ 

Vo C Vi C V 2 C ... C L 

where V, denotes level- j subspace spanned by scaling function <Pj, k W . (*) = 

oi (2V - k) The difference between basis in subspace Vj and V j+ i is the 
wavelet basis'^ (t). In other words, wavelet function ^ t) spans subspace 
w. = Vj+i - V, . It is required that {<p jtk (t) , ^ (*)) - °> for a a PP ro P n 
j, fc, and l. Consequently, wavelet expansion yields 


; (t) = 52 c °- fc ^°- fc W + EE ip j t k (^) 

fc fc i=° 


(33) 
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where 


ij* = (x(t),^w> = / 'difew 


dt 

dt 


In discrete time signal applications, the digital wavelet ran*™ ^(DWT 
operation can be done by using filter banks. The confignratron of DWT filter 
bank is composed of a number of high pass filters, low pass filters, ami dec- 
imators Conversely, to perform inverse digital wavelet transform (IDWT), 
some corresponding high pass filters, low pass filters, and interpolators are 

needed^rder ^ transform a signa i f rom time domain into wavelet domain a 
transformation matrix, T is required. Basically, the dements inthe transfor- 
mation matrix are obtained from the coefficients of high pass filters _and T 
pass filters. The details of constructing the transformation matrix, T can be 
found in [2] and [13]. This matrix is given by 

X (n) = Tx (n) ( 34 ) 

where X(n) is the transformed signal in the wavelet domain and x(n) is a 
sienal in the time domain. 

Now consider 2-band two-dimensional wavelet transform, the same man- 
ner as one-dimension is used to decompose data into subbands except that 
two scaling functions, *>,*(«) M*)' P e ( rform wuvdet transform in the 

horizon and vertical directions. The transformation matrix, T, for two- 
dimensional wavelet transform can be calculated by using those scaling fui 
tions. The transformation equation can also be written as, 


Xll 

XtfH 


T Lt x 

Tl// x 

T//l x 

Thh x 


(35) 

(36) 

(37) 

(38) 


where X p q are images in wavelet domain resulted from performing P opera 
tor in horizontal and Q operator in vertical direction. H and L denote ig 
pass and low pass filtering respectively. And x is an image in time domai . 
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4.1 Wavelet Transform Based Adaptive Filtering 

Wavelet transform is a unitary transform that decomposes a signal by using 
wavelet basis functions. Some work on wavelet based adaptive filtering [1], 
[21 [3] [4], [5], indicated advantages in the wavelet domain over tmcdomm 
The impact of these results to the EDS algorithm comes from subband de- 
properties of the inpot images. Since each decomposed ^ 
contains different frequency components, they are orthogonal to each o • 
“totality oUhe transformation results in faster convergence of the 

EDS - . m 
From the cost function in equation UJ, 

j ( w ) = min {w T Q (n) w - 2w T r (n) + o 1 (n)} 

71 V } „,rp!V L 


we have, 


Q (n) = E [x (n) x (n) T ] 
r (n) = E [d (n) x (n)] 


In the wavelet domain, we transform the input s.gntd x, t^mg ^uat'°n (3^ 
for one-dimensional and equations (35)-(38) for twod.mens.onal. The corre 
sponding correlation matrix Q„ and cross correlation vector r„ are, 


O ... (r).\ = E 


X (r>ML(n) T 


= E [Tx (n) x (n) T r ] 

= TE jx (n) x (n) T j V T 
= TQT t , ( 39 ) 


r w (n) = E [d (n) X (n)] 

= E [d (n) Tx (n)] 

= TE [d ( n ) x (n)] 

= Tr (n) 

The Wiener optimal solution in the wavelet domain is 

w id OP t = Vw 

= [TQT r ] _1 Tr. 


(40) 


17 


(41) 

(42) 


4.2 Image restoration in Wavelet Domain 

In our previous study [8], we have shown the outstanding SNR improve- 
ment in image restoration application using FEDS algorithm in time domain. 
The corrupted noise in testing images is removed by linear predictor using 
( K _ i) x (L - 1) FIR filter. The output of the filter can be obtained as 


K- 1 L - 1 


y ( 711 , 712 ) = S WijX n2-j ) + W °° 

i= 0 j= 0 

where (Lj) + (0,0) , and *(i>i,nj) is the noisy image. Because of the fact 
that images have noil zero mean, coefficient m m is added into the above 
equation. Prom previous section, we know that wavelet transformation can 
be done by using equations (35)-(38). For M-band wavelet transform where 
M = 2, we transform two dimensional data in each level. As a resu. , ere 
are 4 transformed images which contain different information. Each image is 
restored in the wavelet domain separately as shown in Figure 1. As a result, 
we have 4 restored images, Y’ s, in wavelet domain as follows: 


Yll { n li n 2) 

K- 1 L - 1 

= w H Xll ( ni - i5 + w °° 

i = 0 j = 0 
K - 1 L-l 

(43) 

Y lh {ni,n 2 ) 

= yi WjjXLH (7li-t, 7t2-j) + w 00 

t= 0 j=0 
K-X L-l 

(44) 

Yhl ( 711 , 712 ) 

= WijX H L (Til -h n 2-j) + w 00 

i=0 j=0 
K - 1 L-l 

(45) 

Yhh (rii, 7 x 2 ) 

= T E WijXnH { n l-ii n 2-j) + ^00 

t= 0 j=0 

(46) 


The final image in time domain can be obtained by inverse wavelet trans- 
forming of these 4 restored images, i.e. Y L l , Ylh, Y H l,Y H h- 
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Figure 1: Adaptive predictor configuration 

5 Simulation Results 

5.1 Image Restoration 

To demonstrate the performance of the EDS algorithm in the wavelet domain 
we conducted simulations for image restoration in both time and wavelet 
domain. In this simulation, 1-level wavelet transform, with Daubech,e-8 
basis, is applied to lena and clown images. The restoration process using 
FEDS algorithm is performed in each subband image and then e res or 
image is obtained by inverse wavelet transform of those subband 1 »*°red 
images. The simulation results is shown as the improvement of SNR (SN ) 

in Table 3. 

SNR = - 1.8dB 1 SNR = -7 dB ~ 

-Tirnn TWrain I Wave let Domain Time Pom ^hT Wavelet Domam_ 

Lena 9-4 9.9 IL£ 

g aturn 7J 8.7 I 11-0 I J jlL 1 

| Tab | e 3 . Comparison of SNR1 using FEDS algorithm between in lime and wavelet domain. 
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It is obvious that EDS algorithm and FEDS algorithm in wavelet domain 
exhibit improvement over their time domain counterparts. 


5.2 Image Restoration with Lost Samples 

In real world applications in communication systems, received images at the 
receiving end are not only corrupted by noise but also lose some data o 
samples. Therefore, it is necessary to recover those lost data P°“ s P" 
to restoring the images. The configuration of the system is shown Fig . 
x s is an image with lost sample, x' is an image with reconstruction of lost 

samples, and x is a final restored image. , >hed 

To reconstruct the lost samples or data in an image, m 
in fill is used. The reconstruction subsystem is illustrated in lg ■ * _ 

a lost sample image which is an input to the reconstruct, on systenn x is 
obtained by low pass filtering The input image is also, parallely fed to 
hard limiter and rectifier to produce xy Therefore, it .s obviousb - seen that 
x represents sequence of impulses where the lost samples are ocated. x p 
is a low pass filtered version of x p . Division of x S(p and x P(p yie s an iraa 
Withom Lt of samples, x'. Detailed derivation of the reconstruction of lost 
samples in one dimension data can be found in [111. Now we can perform 
noise cancellation of the reconstructed image by passing x to the adaptive 
predictor which is operating in wavelet domain to finally obtain an estima e 



SNR = -1.8dB 

SNR = -7dB 

Lena 

12T 1 

12.5 

Saturn 

8.6 

11.0 


Table 1 SNR I ™„g FEfisIiion.hm b. <R—i« witl. 20% l.H of v,„ P l» ol in, eg,.. 

Table 4 and Figures 4 and 5 are the simulation results when 20% of samples 
in the noisy images are lost. FVom the images and the SNRI numbers, it is 
clear that the developed system achieves significant improvement in image 

restoration. 


6 Multiplierless Filters 


Multiplierless digital filters are important due to their high computational 
speed P In addition, when these filters are implemented in hardware, the elim- 
ination of the multipliers results in a reduction in cost. In two dimensional 
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Figure 2: Image restoration block diagram 



x (n n ) 

p 1 2 


Figure 3: Reconstruction of lost sample subsystem 

(2-D) signal processing, these are very significant advantages because of the 
large amounts of filtering required in such applications as image an ^ 

processing. Several techniques have been developed for desigrn g 
plierless FIR and HR filters during the past decade [18]-[21]. Design me o s 
for 2-D state space multiplierless filters based on a generic aigonthm have 

been given in [18], [19]. Design of multiplierless FIR filters using a McC 

Deen ° . , i m fool foil For the one-dimensional case, 

lan transformation were presented in [20], [21]. tor tne one m 

there are also various techniques proposed for des.gnmg multipl erless i filters 

suS as a design which considers simplified parallel .mplementat.ons of FIR 

filters for signed powers-of-two implementation [ 2 2 ],[ 23 ) tm a ^gi^o mu 

tiplierless IIR elliptic filters based on sensitivity 

erless filters have been designed using periodically shift- variant (PS ) 

[251 The use of periodic filters increases the total number of coefficien s, 
thereby providing the necessary degree of freedom for obtaining power-of- 
two coefficients. If the coefficients are all power-of-two, then they can b 

implemented using simple shifting operations only. , . 2 D 

In this project, we derive the mathematical tools required to design 2 D 
PSV filters with power-of-two coefficients. We start with a given filter m 
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Figure 4: Image Restoration with Lena 


22 



<•) Original Imaga 



(c) Woiay imaga 20* Lo ** Sa»T<p*a 



(b)Nolay Imaga . SNR ■ 7d8 



Figure 5: Image Restoration with Saturn 
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the state-space domain and convert it to an equivalent PSV 
coefficients of this filter can be designed to be power-of-two. An example 
given to illustrate the design methodology. 

6.1 2-D State-Space System 

The state-space model used in this paper is the second model of Fornasini- 
Marchesini (FM) [26], For a linear shift variant (LSV) system, this mode 

can be written as 


x(h + l,k+l) = A l {h,k)x{h,k + 1)+ 
A 2 (h, k)x(h + 1, k)+ 
Ao(h, k)x(h, fc)+ 
B(h, k)u(h, k) 
y(h, k) = C(h, k)x(h, k) 


(47) 


where x{h, k) is an L x 1 state vector, u(h, k) is a scalar mput, y(h k) * * 
scalar output, A 0 (h, k), A 1 (h , k), A 2 (h , k) are L x L state matrices, B(h, k) 

an L x 1 vector, and C(h, k) is a 1 x L vector. 

For the special case when the state matrices are related as 


Ao(h, k) = —A\(h,k)A 2 {h,k) 


and 


A\{h,k)A 2 (h,k) — A 2 {h,k)Ai(h,k), 

the system is called the Attasi model. A filter is called a periodically shift 
varying (PSV) filter with period (P.Q) if all the coefficients o the state 
matrices are periodic with a period (P,Q) he. A(h + mP,k + nQ) - A(h,k) 

for all integers (m, n). ~ TVi^n 

For simplicity, we assume that C(h,k) is a constan vec • > 

the impulse response at position (m,n) due to an input at posi ion { ,j) 
(h(m, n; », j))can be derived in the following way. = 

Let u(i,j) = 6{i,j) , Find y(m,n) for m = i,. + M + 2 ’, n 
i i + 1 7 + 2 ... in a recursive way and then use an induction method to get 
a closed form for an impulse response assuming a zero initial con ltion (i.e. 
x(h,k) = 0 for h <1 or k <l). 
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For m = i, 

(m,n)= y(m,n ) = Cx{i,j ) = 0; 

(m,n)= (ij + 1); y(m,n) = Cx{i,j + 1) = 0; 
(m,n)= {i,j + 2); y(m,n ) = Cx(i,j + 2) = 0; 
(m,n)= (*, j + 3); y{m,n) = Cx{i,j + 3) = 0; 
(m,n)= (i, j +4); y{m,n ) = Cx(i,j + 4) = 0 
For in = i+1, 


(m,n)= (i + l,j) gives 

y(m,n) = Cx(i + l,j) = 0. 

(m,n)= (* + l,j + l) gives 
y(m,n ) = Cx(i + l,j+ 1) 

= CB(i,j)u(i,j) 

= CB(i,j)S{i,j) = CB(i,j). 

(m,n)= (i + l,j + 2) gives 
y(m,n) = Cx(i + 1, j + 2) 

= C\A 0 (i,j + 1 )x{i,j + 1)+ 

A x {i, j + 1 )x(i,j + 2)+ 

A 2 (i,j + l) x (^ + 1) j + 1)] 

= CA 2 (i,j + l)B(i,j). 

(m,n)= (* + 1, j + 3) gives 
y(m, n) — Cx(i + 1, j + 3) 

= C[Ao(i,j + 2 )x(i,j + 2)+ 

A x (i,j + 2)x(i,j + 3)4- 
A 2 (i, j + 2 )x{i + 1) j + 2)] 

= C\A 2 (iJ + 2)A 2 {iJ + 1 )}B(i,j)- 


For m = i+2, 


(m,n)= (i + 2,j) gives 
y(m, n) = Cx(i + 2, j ) = 0. 

(m,n)= (i + 2 ,j + 1) gives 
t/(m, n) = Cx(i + 2, j + 1) 

= C[.<4o(i + 1) jM* + 1) j)+ 

Ai(i -I- 1 ,j)x{i + 1 ,j + 1)+ 

A 2 (i + 1 "i" 2i j)] 

= Cj4i(i 4- \,j)B(i,j)- 
(m,n)= (i + 2, j + 2) gives 
y(m, n) = Cx{i + 2, j + 2) 

= C[A 0 (i + 1, j + 1)®(* + 1> 3 + !)+ 
Ai(i + 1, j + l)x(i + 1, j + 2)+ 
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Ai(i + lJ + l)x{i + 2,j + l)] 

= C[Ao(i + 1 ,j+ 1)+ 

A\{i + 1, j + l)A 2 (i,j + 1)+ 

A 2 (i 4- 1) j 4" l)-4i(i 4" 1) j)\B{i, j)- 
(m,n)= (* + 2, j + 3) gives 
y(m , n) = Cx(i 4-2 ,j + 3) 

= C[Ao(i + 1, j + 2)x(i + 1, j + 2)+ 

A l (i + l,j + 2)x(i + l,j + 3)+ 

A 2 (i 4- 1 ,j 4- 2 )x(i + 2 ,j 4- 2)] 

= C[Ao(i + 1, j + 2)A 2 (i,j + 1)4- 
.4i(i + 1 ,j + 2)A 2 (i,j 4- 2 )A 2 (i,j + 1)4- 
A 2 (i 4- 1, j + 2){A Q {i 4- 1, j + 1)4- 
Ai{i + l,j + l)A 2 {i,j + l)+ 

A 2 (i 4- 1 ,j 4- B)A\(i 4- 1) j)}]^(*? j)- 
For m = i+3, 

(m,n)= (i 4- 3, j) gives 
y(m, n) — Cx(i + 3, j) = 0. 

(m,n)= (i + 3, j + 1) gives 
y(m, n ) = Cx(i + 3, j + 1) 

= C[A 0 {i + 2,j)x{i + 2 ) j)+ 

Ai{i + 2 ,j)x(i 4-2 J + 1)4- 
A 2 (i 4- 2 ,j)x(i + 3 t j)] 

= C[A l (i + 2,j)A 1 {i + l,j)}B(i,j). 

For m = i+4, 

(m,n)= (i + 4, j) gives 
y(m, n ) = Cx(i 4- 4, j ) = 0. 

(m,n)= (i 4-4, j 4- 1) gives 
y(m, n) = Cx(i 4- 4, j 4- 1) 

= C[A 0 (i + 3 ,j)x(i 4- 3 ,j)+ 

A\(i 4- 3 ,j)x{i 4- 3, j 4- 1)4- 
A 2 (i 4- 3 ,j)x(i 4- 4, j)] 

= C[A l {i + 3 t j)A l (i + 2,j) 

A\(i + 1, j)]B(i, j)- 

Since the input is an impulse the output y(m,n) is the impulse 

response h(m, n; i, j) . By induction, we can see that for m > i 4 1 , n > j-\- 1, 

h(m,n\i,j) = CA(m - 1 ,n- 1 )B(iJ) (48) 
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where A(i,j) = In ( N*N Identity matrix ) 

A(m, n) = 0, m < i and n < j 
and for m>i,n>j 

A{m,n) = Ao(m, n)A(m — 1, n — 1) + 

Ai(m,n)A(m — l,n) + 

A 2 {m , n)A(m , n - 1). ( 49 ) 

The above impulse response is for a general LSV system without peri- 
odicity. Now we specialize it to a linear shift invariant (LSIV) system, the 
model of which can be written as 

x(h + 1, fc + 1) = Afx(h, k + 1)+ 

A 2 x{h T 1, fc)T 
A 0 x(h , k)+ 

Bu(h, k) 

y(h, k) = ~Cx(h,k ) ( 50 ) 


with Aq = A\A 2 and AiA 2 = A 2 A v 

The impulse response h(n u n 2 ) of this LSIV system can be derived us- 
ing (48) and (49) with all constant coefficients. The details are omitted for 
brevity. The final result is 

h{ni,n 2 ) = CAy ^ A 2 ^ Bu(ni — l,n 2 — !)• (51) 


6.2 Multiplier- free PSV structure 

It will be shown in this section that a 2-D system consisting of a Hold (which 
increases the signal rate by a factor of (1,Q) ), a PSV system, and a deci- 
mation (which decreases the signal rate by a factor of (1,Q)) connected in 
cascade is a shift- invariant system. Thus we can realize a shift-invariant 
system using a PSV structure as shown in Fig.l. For a multiplier-free 


27 



Figure 6: LSIV equivalent system 


realization, the elements of matrices are restricted in a power-of-two set 
S n = {±2- n ,±2- n+ \...,±l>0}. 

As a first step towards deriving an impulse response of the system ot 
Fig.l, we find the impulse response of the PSV system. . M 

For simplicity, we make the following assumptions: (a) In (47), let A^ti, ) 
be a constant matrix Ai] (b) Ai is commutative with A 2 (h,k)\ (c) A 2 (h,k) 
are commutative among themselves and periodic with period (l,Qh 

Using h(m,n;i,j) in section 2 and replace matrices A lt A 2 and Ac .with a 
constant matrix A u periodic(l,Q) A 2 (h,k) (e.g. A 2 {h + i,k + nQ) -^2 (MO 
for all integer (i, n)) and A 0 (h t k) = -AiA 2 (h, fc) respectively, then the PSV 
impulse response at any position (m,n) due to an input at (i,j) can be denve 


as: 


h(m,n]i,j) — CX(m,n)B(i,j) 

where X (in, n) is given by the following. 

For m = i+1, 

X{i + l,j + l) = lN 

X(i + l,j + 2) = A 2 {i,j + 1) 

X(i + 1, j + 3) = A 2 (i,j + 1 )A 2 (i,j + 2) 

X(i + l,j + 4) 

= A 2 (i,j + l)A 2 {i,j + 2 )A 2 (i,j + 3) 

For m = i+2, 

X(i + 2,j + l) = A 1 

X{i + 2, j + 2) = A\A 2 (i,j + 1) 

X(i + 2,j + 3) 

= AiA 2 (i, j + 1 )A 2 (i,j + 2) 

X(i + 2, j + 4) 

= AiA 2 {i,j + l)A 2 {i,j + 2)A 2 {i,j + 3). 

For m =i+3, 
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X(i + 3, j + 1) = Aj 
X{i + 3,j + 2) 

= A 2 A 2 (i,j + 1) 

X (i + 3, j + 3) 

= A 2 1 A 2 (i,j + 1 )A 2 (i,j + 2 ) 

X(i + 3,j+4) 

= j + l)A 2 {i,j + 2)Mhj + 3 )- 

For m - i+4, 

X(* + 4,j + l) = A\ 

X{i + 4,j + 2) 

= A\A 2 {i,j + 1) 

X(f + 4,j + 3) 

= A\A 2 {i,j + l)A 2 (i,3 + 2 ) 

X(i + 4,j+4) 

= A\A 2 (i,j + l)A 2 (i,j + 2)A 2 (i,j + 3). 

From above, we can see that X(m,n) can be written as 
for Tfi > i + 1 arid n > j + 1, 

n— 1 

X(m,n) = 4T (m) II 

a= j + 1 

Thus, /or m>i + !,«>/ + ! 


h(m,n : i, j) = C+ ( } x 

n— 1 

A 2 {i,a)B{i,j) 

a— j+1 


(52) 


n— 1 

with n +(b«) = /n for n = j + 1. 

Next^we+find the overall impulse response of PSV system of Fig.l. 
With u{n u n 2 ) = 6{n 1 ,n 2 ), we have 

0 Q- 1 . . x 

v(m,n) = d{rn — i.,n — j)- 

i= o j — o u ■ n 

Since the summation over i is from 0 to 0, then l — 

Using the superposition principle, the output can be written as 
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( 53 ) 


y(ni,n 2 ) = w(n u n 2 Q ) 

Q-l 

= y^fc(ni,n 2 Q : 0, j) 

j = 0 

where h{n u n 2 Q : 0, j),0 < j < Q ~ 1 can be ^pressed by using (52) m 
the following for n x > l,n 2 > 1, 

h(n\,n 2 Q : 0,j) =CA^' 1 x 

n 2 Q-l 

J} A 2 (0,a)B(0,j) (54) 

a = j +1 


Q-l . ^ 1 

where J~I ^(O, cl) ^ In when j V 

Now, (54) J can also be written in the following form by using the comma 
tative property of matrices A x and A 2 . For n x > 1 ,n 2 > 1, 


/j(ui , n-zQ 


0,j) = CAl'- l x 

\ri 2 ~l)Q+j 

^ 2 ( 0 , a)) J x 

a = j + 1 
n 2 Q-l 

fl /MO, a) | B(0J) 

t a = ( 1 x 2 — 1)Q+ j +1 


/i(m, n 2 g : 0 ,j) = CA? X 

^(n 2 — 1)Q 

A 2 (0, a) J x 


a — 1 

<3-1 


Y[ A 2 (0, a) j B(0,j) 

\a — j +1 


(55) 
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Notice that the product over one period remains unchanged regardless of 
the starting index. Since the input is an impulse, the output y(ni, n 2 ) in (53) 
is the impulse response of the overall system. Combining equations (53) and 
(55), the impulse response of the PSV system is the following. 

For rii > 1, n 2 > 1, 


h(n u n 2 ) 



Q-l Q-i 

n A 2 {0,a)B(0,j) 

j=0a = j+l 


(56) 


Thus it is shown that the overall system is linear shift-invariant. Note 
that the above impulse response resembles that of the LSIV system in (51). 
Comparing (56) and (51), we get C = C, A\ = A\, 


Q 

t 2 = 

0=1 

Q-l ( Q-l \ 

and B = A 2 (0,a) | 5(0, j). (58) 

j = 0 \a = j+l / 

If we restrict the elements of C and A x to be power-of-two, then we 
have a system whose characterization can be represented by matrices of only 
powers- of- two coefficients. 

There are several ways to find the elements of matrices A 2 and B which 
give the best approximation of the desired LSIV system such as using the idea 
from [25] , or using minimum mean square error criteria. The optimization 
method is under study at this time. 


6.3 Design 

It is assumed that matrices A\ and .4 2 are circulant matrices [25], [30]. That is 
A\ and .4 2 are commutative and all .4 2 are commutative among themselves. 
In addition, we can simplify the problem using the properties of circulant 
matrices as shown below. 
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Any circulant matrix A of order N x N can be diagonalized as A = F*AF , 
where A is the diagonal matrix of eigenvalues, F is a DFT matrix whose k\ th 
element is exp , 1 < k,l < N and F* is the complex 

conjugate of F which has property that F* = F _1 . We diagonalize both 
sides of (57) to get 

0 

T 2 = Y[M0,a) (59) 

a = 1 

Let FB = P and FB = B. Then from (58) we have 

0-1 / 0-i \ 

P = J2 F ( U M0,a)\F*FB(0,j) 

j = 0 \a = j+ 1 / 

0-1 / Q-i \ 

= E II A 2 (0,a))B(0,j) (60) 

j = 0 \a = j +1 / 

Using the criteria based on a square error in time domain, 

{h d (n u n 2 ) - h{n u n 2 )) 2 (61) 

ni ri2 

where (rn,n 2 ) € the region of support of h d (ni, n 2 ), which is a desired 
impulse response. 

Assume that a desired filter is given in a state-space form as in(51). 
for ni > 1, n 2 > 1, 

h d (n u ni) -CA^'W^B (62) 

otherwise, h d (rii,n 2 )=0 _ _ 

Using (56) for h(n\,n 2 ) with C = C, A\ = A\, and since all matrices 
are circulant, we can diagonalize (62) and (56) using a DFT matrix. Let 
~CF* = C and combine with (59) and (60) to get 
h d {m,n 2 ) = ~CF*Fl[ (ni ~ l) F*FT 2 {n2 ' l) F*FB _ 

= ~CF\FT^ ni ~ l) F*){FA} n2 ~ l) F*)FB 

h{n u n 2 ) = ~CF*FT 1 (ni ~ 1) F*F 
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( n^(o> a )) Ba " 1F ‘ F 

a — l 

Z ( rj 1 A 2 (0,a))F*FF(0,j) 

j — 0 a = j +1 

= ~CF*(FT l {ni ^ 1) F*) 

(UFA 2 (0,a)F*) n2 ~ 1 

a= 1 

Q E( Q n F^2(0,a)F*)FB(0,j) 

j = 0 a = j +1 

= CM (ni-1) (n A 2 (0,a)r- 1 

a = 1 

Q E( JI A 2 (0,a))F(0,i) 

Thus the square error in time domain (61) can be written as 

I \ 

(flA 2 ( 0 ,a)r-‘ 

& = zL Z_> a = 1 

ni ri2 Q — 1 Q— 1 ~ 

£ ( n A 2 (0,a))B(0J)) j 

Note that instead "of finding combination of circulant matrices A 2 (0, a), 
we find only combinations of diagonal matrices of eigenvalues of A 2 (0,a). 


6.4 Example 

Consider a LSIV 2-D Attasi’s model filter using circulant matrices as follow- 

0.5 -0.5 0.125 -0.125 

-0.125 0.5 -0.5 0.125 

0.125 -0.125 0.5 -0.5 

—0.5 0.125 -0.125 0.5 

0.5 0 -0.015 0.25 

_ 0.25 0.5 0 -0.015 

A2 = -0.015 0.25 0.5 0 

0 -0.015 0.25 0.5 

B = [ 1 0.39 -1 0.45 ] T 
C= [ 1 -i -i 1 ] 
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Using the criteria (61) above with periodicity Q - 2 and elements of 
, . c _ / 1 zd 0 - - li* the designed matrices 

matrices belong to set b 2 - \ -1 > 2 > 4 ’ u ’4>2’ /’ 

are obtained in the following: 

A\ = Ai and C = C 

” 0.5 0.5 


A 2 (0,0) = 


^2(6) 1) 


-0.5 

0 

0.5 

0.5 

0.5 

0.5 

0 


0.5 

-0.5 

0 


0 

0.5 

0.5 

0.5 


B ( 0 , 0)=[-1 -1 

B{ 0,1) = [1 "0.25 


0.5 

0 

0.5 
0.5 

-0.25 


0 

0.5 

0.5 

-0.5 

0.5 

0.5 

0 

0.5 


-0.5 

0 

0.5 

0.5 


-ir 
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This design yields an error E = 0.0108. The frequency responses of the de- 
sired LSIV filter and the designed multiplierless filter are quite we ma c le 

as shown in Fig. 7. 





Fig. 7a.) Desired magnitude response. 
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Fig. 7b.) Designed magnitude response. 

T Conclusion 

In this project, we have developed a new adaptive filtering algorithm for im- 
age restoration and compression. The algorithm is called the Fast Euclidean 
Direction Search (FEDS) method. We have presented the derivation of this 
algorithm and all its virtues. This algorithm has then been implemented in 
the wavelet domain in conjunction with an algorithm for lost sample recovery. 
The overall system is then shown to perform efficiently for image res ora ion. 
We have also presented some new results on the design of 2-D mu tip ler ess 
filters. These filters will be used in the future for image and video processing 

applications. 
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